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Abstract 

In order to get an estimate of the homogeneity of the distribution of matter in a fast 
hadron or nucleus, we compute the correlations of the saturation scales Qa between 
different points in impact-parameter space, in some specific saturation models. We 
find that these correlations are quite strong: The saturation scale is nearly uniform 
in domains whose sizes scale like exp [const x ln^(l/a^)] /Qs, which means that the 
density of gluons should not fluctuate significantly over regions of that typical size. 
We expect these conclusions as well as the explicit analytical expressions we obtain 
for the correlations to be true also for full QCD in appropriate limits. 



1 Introduction 



In the high-energy regime of QCD, an interesting new phenomenon is expected 
to show up: parton saturation [l;2j. Saturation changes qualitatively the usual 
equation for the evolution of scattering cross sections with the energy of the 
reactions, namely the so-called Balitsky-Fadin-Kuraev-Lipatov (BFKL) equa- 
tion [3J , by introducing nonlinearities. In turn, due to these nonlinearities, the 
intrinsic stochasticity of partonic evolution may start to have a sizable effect 
on observables. 



The basic nonlinear evolution equation beyond the linear BFKL equation is 
the Balitsky-Kovchegov (BK) equation [1|5] (or, alternatively, the Balitsky- 
Jalilian Marian-Iancu-McLerran-Weigert-Leonidov-Kovner (B-JIMWLK) equa- 
tions [4|6j) which however neglects the stochastic effects alluded to before, and 
is therefore a kind of mean-field approximation. While a formulation which 
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would include all stochastic effects has not been fully achieved yet (the most 
recent advances may be found in Ref. [7]), it is believed that the complete evo- 
lution has a lot in common with some react ion- diffusion processes described 
by equations of the Fisher-Kolmogorov-Petrovsky-Piscounov (FKPP) type [S]. 
At the mean-field level, this analogy is a formal identity between the FKPP 
equation and the BK equation in the so-called diffusive limit and assuming 
uniformity in the transverse space (i.e. the gluon distribution is assumed to 
evolve in the same way at all points of the impact-parameter space) |9]. Be- 
yond the mean-field approximation, the conjecture made so far for the full 
problem is that at any fixed position in transverse space, the rapidity evolu- 
tion of say the gluon content of a hadronic object is like the time evolution of 
a one-dimensional reaction-diffusion process, whose space variable would be 
the logarithm of the transverse size (or momentum) of the gluons. Generally 
speaking, the dynamics of such systems is described by equations equivalent 
to stochastic extensions to the FKPP equation. (The first ideas on how to 
go beyond the BK equation were presented in Ref. [10]; the equivalence with 
reaction-diffusion processes was conjectured in Ref. [11]; deeper insight can be 
found in e.g. Ref. [12], and a review in Ref. [13].) 

While this analogy is useful to find asymptotic properties of observables that 
involve one unique impact parameter, so far little is known about the corre- 
lations and the fluctuations of the gluon distribution between different points 
in transverse space, which would show up in observables that probe several 
points in impact-parameter space simultaneously. 

The correlations of the gluon number densities at different points in transverse 
space were computed exactly in Ref. [H] in the context of the systematic ap- 
proximation to full QCD provided by the color dipole model flS], which is an 
accurate representation of the physics described by the BFKL equation but 
which does not take into account saturation effects. In Ref. [IH], a calculation 
was done in a theory with full saturation based on a similarity with the Liou- 
ville gravity. But that calculation was valid only up to distances of the order 
of the inverse saturation momentum. On the numerical side, on one hand, 
the BK equation was solved taking into account the full impact-parameter 
dependence [17], and on the other hand, the relevance of the one-dimensional 
stochastic FKPP equation at each fixed impact parameter was tested in toy 
models for QCD evolution beyond the mean- field approximation [18j. 

In this paper, we would like to investigate how the saturation scale varies in 
impact-parameter space in the presence of both saturation and fluctuations. 
Our method will consist in proposing a simple toy model which contains the 
main physical features of QCD, which may be implemented as a Monte-Carlo 
event generator and for which analytical calculations will be possible. In these 
respects, our approach follows the one developed in Ref. [iSj, but while the 
latter work was purely numerical, our main results will consist in analytical 
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expressions of the correlation of the saturation scale between two points in 
impact-parameter space, as a function of the distance between the points and 
function of the rapidity. 



The model is introduced in the next section. We then provide the derivation of 
the analytical expression for the correlations. Finally, we check our calculations 
against numerical simulations of different versions of the model. 



2 Toy model 

2. 1 Construction 

The model that we introduce here is a simplified version of the model proposed 
in Ref. [18]. 

The starting point is the QCD color dipole model [IS], supplemented with 
some ad hoc saturation mechanism which limits the number of dipoles in 
any given phase space cell. The dipole model accurately represents the QCD 
evolution in the high-energy regime and in the limit of a large number of colors. 
It provides an equation for the change of the density of gluons (represented 
by a set of color dipoles of different sizes and positions in the two-dimensional 
plane transverse to the flight axis of the hadron) inside a hadron, when rapidity 
is increased. The basic process from which the evolution is built is the splitting 
of a dipole represented by its two endpoints (xq, xi) into two dipoles (xq, X2) 
and (xi,X2) with the rate [T3] 



As usually, a = asNc/n, where is the strong coupling constant, A^^^ the 
number of colors and Y is the rapidity. The splitting of dipoles is a linear 
process, which generates the BFKL equation when averages over dipole con- 
figurations ( "events" in an experimental language) are taken. The rapidity Y 
is an effective evolution time. 

When the rapidity becomes very high, then gluons and thus dipoles may start 
to interact among themselves, which induces nonlinear terms in the evolution 
equations. The effect of these interactions is to tame the growth of the phase- 
space number density of dipoles as soon as it reaches ~ l/'^s) which would 
otherwise be exponential with the rapidity. The precise mechanism for these 
effects is still not known in QCD, but the main observables should be quite 
independent of these details. 



d{aY) 



dP 
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With respect to QCD, we assume the following simplifications: (i) Dipoles 
evolve by giving birth to one dipole of half size (the left or the right half of 
the parent dipole), or to one dipole of double size (in such a way that the 
parent be the left or right half of its offspring) at some fixed rates, (ii) dipoles 
do not disappear in the evolution, that is to say, the parent dipoles are not 
removed, (in) the positions and dipole sizes are discrete, and (iv) the configu- 
ration space of the dipoles is a line instead of the full two-dimensional space. 
We thus give up two main properties of the QCD dipole model: The collinear 
singularities, which cause the dipole endpoints to emit an arbitrary number 
of dipoles of arbitrarily small sizes, and the continuous and two-dimensional 
nature of the dipole sizes and positions. The first simplification is the diffusion 
approximation, which has been studied in the context of BFKL physics (see 
e.g. Ref. [19j), but which was not assumed in Ref. ^J. The second simplifi- 
cation was instead already assumed in [TS]. These model simplifications may 
introduce some artefacts, but that we believe are under control, and many 
results which we will obtain within such simple models are likely to apply to 
QCD since they will not depend on the details. 

Let us now specify completely the model. According to the evolution rules 
given above, starting from a dipole of size 1, the sizes of all dipoles present in 
the system after evolution are powers of 2. In practice, we shall only consider 
fractions of 1, i.e. the sizes may be written as 2^"^^, where k > 1. For each 
value of fc, there are 2^~^ possible values of the position h of the center of 
the dipoles: h = -\ + 2-^ -| + 3 x 2-^ ■ ■ ■ , | - 3 x 2-^ \ - 2"^ Let us 
number these bins by the index < j < 2'^"^ — 1 running from the negative to 
the positive positions. The model may be represented as a hierarchy of bins 
that contain a discrete number of dipoles, see Fig. [H Note that to any given 
impact parameter h between — | and | corresponds one unique bin at each 
level of size. For example, at position h = — |, one sees the bins {k = l,j = 0), 
{k = 2,j = 0), {k = 3,j = 0) etc... At position —0.2, one sees the bins 
{k = l,j = 0), {k = 2,j = 0), {k = 3,j = 1) etc... More generally, at position 
— I + Ab, one sees {k, [Ab x 2'^^-'^]), where the square brackets represent the 
integer part. 

During the rapidity (or time) interval dt, a dipole in the bin {k,j) has a 
probability adt to give birth to a dipole in the bin {k + 1, 2j), adt to give a 
dipole in the bin {k + 1, 2j + l), and (3dt/2 to give a dipole in the bin (fc — 1, j'/2) 
if j is even and {k — 1, {j — l)/2) if j is odd. Note that dt may be infinitesimal 
(which is generally speaking convenient for analytical calculations), but also 
finite (which is convenient for numerical simulations). 

As for the saturation mechanism, we assume the simplest one: We veto split- 
tings to bins which already host the number of dipoles. 

We can consider that the number density of "gluons" of a given size seen at 
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Fig. 1. The hierarchical structure of the model. Each box represents a bin which 
may contain up to N dipoles of given sizes (vertical axis) and positions in impact- 
parameter space (horizontal axis). The conventional numbering of the bins that we 
have chosen is also shown for = 1, 2, 3. 



one impact parameter is proportional to the number of dipoles in the corre- 
sponding bin {k,j). As rapidity is increased, the occupation of the bins with 
low values of k gets higher until the number of objects they contain reaches 
N. The subsequent filling of the bins indexed by larger values of k (smaller 
dipole sizes) can be seen as the propagation of traveling wave fronts at each 
impact parameter, with possibly complicated relationships between them. The 
(logarithm of the) saturation scale X{b, t) at impact parameter b is related to 
the position of the front seen there at time t. There are several equivalent 
ways to define the position of the front. It could be, for example, the largest 
value of k for which the number of objects becomes some given fraction of N. 
(Later, we will use a slightly different definition). 



2.2 Basic features of the model 



Let us denote by ?^(fcj)(t) the number of dipoles present in bin {k,j) at time t. 
Then, according to the rules given above, we can write the following stochastic 
evolution equation: 



n{k^j){t + dt) — min 



{t)-r5 



/3/2 



, (2) 
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where the 5^f^ j-^ are drawn according to the binomial distribution 



Proba 



This is a rather comphcated equation which we do not know how to solve 
except numerically. 

This model does not a priori look like a stochastic FKPP model. We may 
assume uniformity in impact parameter: This would amount to imposing the 
same and 5^/^ respectively for all j at any given k. In this case, the model 
would be projected to the FKPP class, but by definition, this would wash 
out the fluctuations between the different impact parameters. This simplified 
model, that we call "FIP" (for "Fixed Impact Parameter" , since effectively, the 
model is completely defined by a single impact parameter) in the terminology 
of Ref. [18J, is nevertheless useful since it provides a benchmark to evaluate 
how the fluctuations between different impact parameters may alter the FKPP 
picture. In this paper, we will rely on (and check again in the case of our model) 
the conclusion reached in Ref. p2j that thanks to saturation, locally at each 
impact parameter, the full model is still well-described by a one-dimensional 
FKPP equation, and the fluctuations between different positions in impact- 
parameter space do not qualitatively change the picture. 

Let us first apply the well-known treatment of FKPP equations to the FIP 
case. We know that the large-rapidity realizations of the model are stochastic 
traveling waves, whose main features can be determined from a simple analysis 
of the linear part of the evolution equation. In this model, only the number 
of dipoles Uk in the bins sayXfc,0) (i.e. at impact parameter — |) is relevant. 
The evolution equation read^lJ 



Uk (t + dt) = min 



iV,nfc(t) + 5,%(t) + CiW 



(4) 



The mean-field (or Balitsky-Kovchegov) approximation to the evolution leads 
to the equation 



Ukit + dt) 



mm 



N^rikit) +adtnk-i(t) + f3dtnk+i(t) 



(5) 



where the nk are now real functions of k. The linearized equation (equivalent 
to the BFKL equation) is simply obtained by discarding the "min" in the 



^ We could also write 26^^^{t) instead of the last term in Eq. @. (This may even 
be a more literal implementation of the FIP approximation). But this would not 
make a large difference, which anyway, we would be unable to capture analytically. 
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previous equation: 



nk{t + dt) = nk{t) + adt nk-i{t) + [5dtnk^i{t). (6) 

From standard arguments, we know that for asymptotically large t and iV, 
the velocity of the wave front, that is the time derivative of the position X{t) 
of the front, is given by piT|2Uf21] 

^^o = ^ = x'(7o), (7) 

where x{l) is the eigenvalue of the kernel of the linearized evolution equa- 
tion ([6]) corresponding to the eigenfunction e"'''^, namely 

x{-i) = j^hi(l + adte' + l3dte-^), (8) 

and 7o minimizes xil)/!- We recall that dt may be finite or infinitesimal, in 
which case Eq. ([H]) is to be understood as the derivative of In. The shape of 
the front is a decreasing exponential to the right of the saturation region, 

nfc(t) ~ Are-^«('=-^(*)), (9) 

extending to A; — )■ oo for large times. 

The corrections due to the discreteness of are known. Taking the latter into 
account, the front has now a finite extension, of the order of 

Lo ^ (10) 

7o 

Indeed, typically, it cannot extend further than the point where n^it) ~ 1, 
and Eq. fllOp then follows from Eq. Q. Its mean velocity reads P^PH] 



^^BD = X (7o) - ^ j2 ■ (11) 



This velocity was obtained in a still deterministic (mean-field) approxima- 
tion with appropriate cutoffs to mimic saturation and discreteness. Taking 
furthermore fluctuations into account, the velocity becomes 



'( \ 7r^x''(7o) i.^. 
v = X (7o) z ^, (12) 

27o(Lo + 4lnLc 



and realization-to-realization fluctuations are characterized by the following 
cumulants of the position of the front: 



[ra— th cumulant] -k"^ x" {lo)n\C,{n 
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where a priori the model leads to L = Lq, but empirically, a better fit to 
the results of the numerical simulations is obtained by adding a subleading 
correction of the form L = Lq + 31nLo/7o + const. In particular, the diffusion 
constant of the front, that is to say the slope of the t-dependence of the 
second-order cumulant reads 

^ = ^- (14) 

These expressions are valid in the limit of large L ~ Lq, i.e. for exponentially 
large values of N. 

There exists a subclass of these models that may be reduced exactly to a 
collection of FKPP models. Let us set /3 = 0, that is to say, authorize only 
splittings to smaller-size dipoles. Consequently, there may not be any influence 
of the bins at any size level k on the content of the bins of level less than k 
(i.e. of larger sizes). Then, at each impact parameter, one has a FIP model, 
i.e. a model of the one-dimensional FKPP type. However, the relationship 
between the different copies of FIP models is not trivial, since part of the 
evolution is common between different impact parameters. Even for /9 7^ 0, we 
believe that this minimal model, consisting in considering two one- dimensional 
systems appropriately correlated, is a good approximation to the full model. 

In our investigations, we will have in mind the latter class of models, and 
we will check numerically that for more general models for which /3 7^ 0, the 
results that we shall obtain are not significantly altered. 



3 Correlations in the hierarchical model 



Our aim is to study the correlations between the point at position b = in 
transverse space (left edge of the system, see Fig. [1]) and the one at position b = 
— I+A6 with < Ab < 1. We calculate the average of the squared difference of 
the positions of the front between these points, which is formally related to the 
two-point correlation function of (the logarithm of) the saturation scales, and 
which we deem a good estimator of the spatial fluctuations of the saturation 
scale. In the hierarchical model, all bins with index k less than or equal to 
kAb = 1 + [— log2 A6] (the notation "[■ ■■]" stands for the integer part) and 
j = overlap both impact parameters, and thus the dipoles of size larger than 
2-kAb seen at these points are exactly the same. For k > k^b instead, the bins 
seen at the two points are distinct and nonoverlaping. So in particular, in our 
model with /3 = 0, as soon as the position of the front at one point or at the 
other is larger than A;a65 that is to say, as soon as there are of the order of 
dipoles in the bin {k^b^j = 0), then the evolutions are completely uncorrected 
at the two points in the corresponding bins. (We expect that for finite /3 of 
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order 1, the discussion would not be qualitatively changed.) This matches to 
the picture that we may infer for the QCD dipole model: The dipoles at two 
positions in impact-parameter space separated by a distance larger than the 
typical saturation scales in that region evolve (almost) independently towards 
larger rapidities. Note that choosing pairs of points around impact parameter 
0, one with positive impact parameter and another one with negative impact 
parameter, would not satisfy this property, due to the rigidity of the sizes 
and positions of the dipoles. Indeed, these two points would decorrelate very 
soon in the evolution since their common ancestors necessarily sit in the bin 
(k = 1, j = 0), see Fig.m 

As a consequence of these features of QCD reproduced in the toy model, 
studying two-point correlations between points in impact-parameter space as 
a function of their distance Ab and of the time (=rapidity) t is equivalent to 
studying the time dependence of the correlations of the saturation scales of 
two realizations of the model whose evolutions are identical until the tip of 
the front reaches k/^^. On the average, it takes a time ^Afe = (^Ab — v 
being the mean velocity of the individual fronts, for the front whose tip is at 
^Ab = 1 at the beginning of the evolution to have its tip at k^b- Then the bins 
such that k > k^b evolve independently between the two realizations over the 
remaining time interval 

At = t- tAb, with tAb = [~ ^"^2 . (15) 

V 

Note that this is very close to assuming that the realizations are identical for 
t < tAb and completely uncorrelated for t > tAb- 

From this discussion, we see that the basic input of our calculation will be the 
mechanism for the propagation of a FKPP front. We will review it in the next 
subsection, then we will proceed to the formulation of the calculation of the 
correlations. 



3. 1 Short review of the mecanism for stochastic front propagation 

In this section, we review the recent progress in the understanding of the 
mecanism for front propagation [22|10|23j . which has eventually led to the 
formulation of a model, and from which Eqs. f lT2|) . f|T3|) follow [23] . 

Instead of solving the nonlinear problem, it was proposed P^fTU] to replace the 
nonlinearities by absorptive boundaries, and to treat the evolution equation 
as a linear (branching diffusion) equation between these boundaries. There are 
two types of nonlinearities: (i) The saturation condition that keeps the number 
of particles in each bin less than or equal to N and (ii) the very discreteness 
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of this number of particles, whose effect mainly shows up in the region where 
n ~ 1. The position of the boundaries is adjusted in such a way that the 
distance between them be the size of the front fllUp . namely Lq = lniV/70, 
and that the large-time solution of the diffusion equation be stationary. This 
procedure leads to the expression of the front velocity for large given in 
Eq. f lTTj) . but it does not predict the fluctuations of the position of the front. 

In order to incorporate the latter, we add to this picture the possibility that 
there be one (or a few) particles randomly sent ahead of the tip of the front |23] , 
at a distance 6. We attribute to this event a probability per unit time 

p{6)d6 = Cie-"">^Q{6)d6, (16) 

in the continuation of the shape of the front (Q solution in the large-t and 
large- limit. (Ci is some constant, undetermined at this stage). In Ref. [23] . 
we treated such fluctuations as localized extra weights of appropriate "mass" 
a unit or so to the left of the right boundary. We then computed the effect 
of this weight on the position of the front at large time, and found a forward 
shift of the position equal to 

i^(,) = _,„(l + C,^j, (17) 

C2 is another constant. At large times after a fluctuation has occurred, the 
front relaxes to its mean-field shape. We assumed that the relevant fluctuations 
are rare enough in such a way that the front has time to completely relax 
between two fluctuations, which turns out to be true for Lq ^ 1. 

From the two quantitative elements (I16p and f lT7|) . together with the solution 
of the stationary diffusion problem between the boundaries flTT]) . we may write 
down an effective theory for stochastic front propagation, with however two 
unknown parameters, namely Ci and C2. While we were not able to determine 
Ci and C2 separately, it is the product C1C2 that appears in all cumulants 
of the position of the front, and in particular in the correction to the velocity 
induced by fluctuations. Assuming that the velocity of the front taking into 
account the fluctuations is the velocity of a mean-field front whose size is 
extended by 31nLo/7o with respect to Lq (i.e. with the substitution Lq — )• 
Lq + 31nLo/7o), we got a determination of the product C1C2: 

CiC2 = n\"{^o). (18) 

The effective theory for front propagation then leads to Eqs. (fT2l) and (fT3|) . 
As we will see below, the only new ingredient that will be needed is the 
time dependence of the shift of the front R, whose large-time asymptotics is 
Eq. (dZ]). 
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3.2 Formulation of the calculation of correlations 



In line with the above discussion, we wish to compute the correlations of the 
position of two fronts whose evolutions are identical for t < ^Afe and uncor- 
related for t > t^b- Note that strictly speaking, we would need to keep the 
content of all bins k < kAt identical between the two realizations at all times, 
even after time tA6- But these two formulations give quantitatively similar 
results. 



Let us introduce X{to,t) the position of the front at time t in the frame in 
which X{to,tQ) = 0. We focus on what happens slightly before the initial 
time to- On one hand, X(to — dto, t) = X(to, t) + v^ndtQ if no fluctuation has 
occurred between times to ~ dtQ and t^, on the other hand, X(to — dto^t) = 
X{tQ, t) +VBDdto+ R{t —to, S) if a fluctuation has occurred at a position 6 ahead 
of the front (which happens with probability p{S)dS dto). It is straightforward 
to write an equation for the generating function of the cumulants of X: 

_ i_ln/e^^(*«'*)) = XvBB + f d6p{6) (e^^(*-*'"^) - l) . (19) 

(JjLq J 

One now considers two such independent fronts and add up the generating 
functions. One gets 

d6 p{6) (e^^(*-*°'^) + e-^^^*-*^''^) - 2) . (20) 



Expanding for A close to 0, the coefficients of the second power of A obey the 
equation 

^ ((Xi - X2Y) = 2 J d5p{5)R\t - to, 6), (21) 

where we have used the fact that Xi and X2 are independent random variables 
for t > to, and we have traded to for t in the derivative, taking advantage of 
the fact that both Xi — X2 and R only depend on t — to- In practice, to will be 
equal to t^b, the time at which the tip of the single front reaches k^b. From 
Eq. ([T5|) . this time is [— log2 A6]/t>. 



We see that the basic ingredient is the time evolution of the shift of the front 
due to a forward fluctuation. This shift was given in Eq. (fT7|) in the limit of 
large times. We have to repeat the steps that led to Eq. ( IT7|) keeping however 
track of the full time dependence. 
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3. 3 Effect of a fluctuation on the position of the front 



The problem amounts to solving a diffusion equation between two fixed ab- 
sorptive boundaries, with various initial conditions. We shall discuss the scaled 
dipole number u{t, k) = nk{t)/N, and for the simplicity of the formalism, con- 
sider that it is a function of a real variable k. 



Let us write the general branching diffusion equation: 

dtuit,k) = xi-dkHt,k), (22) 

where x is an appropriate kernel that encodes the linear evolution of the 
dipoles. (The operators that appear here are denoted as differential operators, 
but they could also be finite differences as in Eq. ([6]). Analytical calculations 
are usually easier with differential operators). In the case of QCD, x{~dk) 
would be the BFKL kernel (that can easily be deduced from Eq. (II])), t ~ aY, 
~ In (where r is the size of the dipoles), and u would be the scattering 
amplitude. Further, we define 70 to be the eigenvalue of x that satisfies x(7o) = 
7ox'(7o)- Following Ref. [23], we write the ansatz 

= (H2^. . (23) 

L is the size of the front, which is essentially Lq = In for large A^. When 
X is expanded to second order around the eigenvalue 70, then ip obeys the 
partial differential equation 

dy^ = \dli^ + ^^^(x'(7o) - X'(t))^, (24) 
where we have defined 

2x"(7o)t , k-X{t) 

y = — — P ^ — L — ■ ^ ^ 

We have only kept the dominant terms for large L. We see that x'(7o) ~ 
X'{t) has to scale like for all terms of this equation to be relevant. The 
coefficient must be chosen in such a way that in the large-?/ limit, there is a 
nontrivial stationary solution. We shall use the already known result ( II ip to 
write 



X'{t) = x'ilo) - ^^^^ + 0(1/ L') (26) 



and check a posteriori that it is the correct expression. Equation I HM then 
becomes 

dy^ = ^^pV + (27) 
up to higher-order terms when L is large. 
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We shall admit that the saturation of the number of dipoles and the stochastic- 
ity may be appropriately implemented, within well-controlled approximations, 
by respectively an absorptive boundary at p = and another one at p = 1 
(which corresponds to a distance L between the boundaries in fc-coordinates, 
i.e. to the natural width of the stationary front which travels at the velocity 
fBo)- The boundary conditions formally read 

^{y, p = 0) = and ij{y, p = 1) = 0. (28) 



Let us now discuss the initial condition. The fluctuations that will generate 
the front-shifts R responsible for the corrections to vbd are dipoles sent at a 
distance 6 ahead of the tip of the steady front. Formally, in order to describe 
such a fluctuation occurring at time t = 0, we would write m(0, /c) = 5{k — 
X{0)-L-6)/N, i.e. fromEq. ([23D, ^p{0,p) = e^°^5{p-l- 5 / L) / L"^ . But since 
there is an absorptive boundary at p = 1, we had better move the fluctuation 
to the left of the right boundary. This should not qualitatively change the 
problem as long as L ^ 1. Thus we write 

V'(y = 0,p)=5(p-l + a) — , (29) 

where a is a constant of order and therefore a ^ 1. 

The solution of Eq. f l7r|) with the conditions and (12^ reads 

ipsiVyp) = (—1)"'*'^ sin Tina sin vrnp e * . (30) 

^ n=l 

The subscript 6 recalls that this solution was obtained starting with an initial 
condition consisting in one dipole at a distance 6 ahead of the deterministic 
front. 

We first discuss the stationary solution. We see that for large y, the higher 
harmonics are suppressed exponentially with respect to the fundamental mode 
n = 1, which gives the following contribution: 

2e7o5o 27rae'^°^° 
i^So (y, P) = -JY~ — J2 — ^P- (3^) 

Thanks to the choice (1261) for X'{t), this solution has no y dependence, and 
leads to a stationary u in the frame of the front. The expression ( 13T|) is inde- 
pendent of the initial condition except for the overall normalization. The value 
of 6q, which characterizes the initial condition, will be adjusted later. Undoing 
the changes of variables which trade u for ip, k for p and t for y (Eq. (^3^), 
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the stationary solution us^ reads 



L sin 



■K{k-X{t)) 



(32) 



We require that uo{t,k) ~ 1 for = X{t) + aL, where aL is a constant of 
order 1. This condition is satisfied if we set 6o ~ 31nL/7o. Indeed, with this 
choice, 

us^it, X{t) + aL) ~ 2'K'^aaL^e-^°''^. (33) 
Since aL and aL are constants, the right-hand side is just a number of order 
1. 



We now add a fluctuation to the stationary front, that is to say, an extra 
particle at a position 5 ahead of the tip of the steady front. The solution of 
the diffusion equation is the superposition of the large-time stationary solution 
given by Eq. (!32|) with 5q = 3 In L/70, and of the solution us of the diffusion 
equation with the generic initial condition characterized by S (see Eq. (130|) ). 
up to a multiplicative constant C2 of order 1 that we do not control in this 
calculation, since it certainly depends on the detailed shape of the fluctuations. 
We write 



u{t, k) = usa{t, k) + C2Us{t, k) 

= e-^«('=-^W)L[^5„(y,p) + C2^,(y,p)] 



(34) 



up to the replacement of the variables by their expressions fl2Sl) . The presence 
of the second term alters the shape of the front (the front eventually relaxes 
back to the sine shape in Eq. (132|) ). see Fig. [21 But of course, we want to keep 
the normalization condition for m, namely for some appropriate value of k, u 
is required to equate Eq. (133!) at all t. This is possible by shifting the value 
of k at which we enforce the normalization condition from k = X{t) + aL to 
say k = X(t) + aL + R{t, 6). This is equivalent to shifting the position of the 
front X{t) X{t) + R{t, 6). Equation then leads to 



u{t, X{t) + aL + R{r, 6)) = 2n^aaL'^e-^'>^^-^'>^^^'^^ 



a 



2n'^aa L 



(35) 



Equating the right-hand sides of Eq. ( j35l) and Eq. ( j33l) . we get 



R(t, 6) = — In 
7o 



2'n'^aaL 



(36) 



where only the lowest orders in a, a in the expansion of ips must be kept. 
With the help of Eq. (130|1 . it is then straightforward to arrive at an explicit 
expression of R. 

Note that by performing the shift X — )■ X + /2, we have actually added a 
t-dependent term to X'{t) in Eq. which was initially assumed to be a 
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Fig. 2. ipSgiuTp) + ips{y,p) [see Eqs. (f3T]) . (j30]l ] for different values of the reduced 
time variable y after a fluctuation of size 5 = 5 has occurred at y = 0. The size of 
the front is L = 10, and a = 0.1. We see how the fluctuation, initially localized at 
the tip of the front, gets smeared uniformly over the width of the front as y gets 
large. Eventually, a small forward shift X ^ X + R would be needed in order to 
absorb it and recover the stationary front. 



constant. However, since the time scale for the variations of R is L^, R' is 
suppressed for large L, and thus Eq. (|271) is not modified within our approxi- 
mations. 

It is interesting to note that ips is related to some Jacobi function [23]. Since 



i^^(z\q) = l + 2^(-l)"cos(2n2)g"'. 



(37) 



n=l 



we may rewrite Eq. (!30ll as 

g7o5 1 



My,p) 



2L2g 



^^4 



'7r(a + p) 



^4 



TT 



(a-p) 



(3^ 



The notation 

q = e ~ = e (39) 



T^x"(7n)t 



has been introduced. Using Eq. (!36|) and performing the appropriate expansion 
for small a and a, we arrive at an expression for R{t, 6) in terms of the ■i?4- 
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function which is particularly compact: 




(40) 



with 



5,^4(0|g) = 2 5:(-l) 



(41) 



n=l 



It is actually quite natural that the Jacobi theta functions appear, since the 
latter are defined as solutions of the one-dimensional heat equation with pe- 
riodic boundary conditions. 

We turn to the analysis of the obtained result. First, for large only the 
fundamental mode contributes to ips, and it is clear that Eq. ( 136|) reduces to 
Eq. f|T7j) . Looking back at Eq. fl30l) . we see that higher harmonics would give 
a series of exponentially decreasing corrections. But at a finite time, a large 
number of modes have to be taken into account, typically all modes such that 



n < (L/vr) Y2/x"(7o)t. A few low-lying modes are not enough to describe the 
small-time behavior. Instead, it is a saddle point (in an appropriate integral 
reformulation) that dominates the sum ( !30|) . In this regime, it would be useful 
to find a way to write the series of harmonics such that at asymptotically large 
only the first term contributes instead of the whole series. This is actually 
possible using the Poisson summation formula 



In order to get i?(t, 5), we need the value of ips aX p = a. Hence we choose 

g70<5 2 

f{x) = —smnxasmiTxaq^'~^e^'^^. (43) 

We then perform the integral over x in the r.h.s. of Eq. ( l42l) . Introducing 
7+ = a -|- a and 7_ = a — a, we get the following expression for tps'- 



Since we eventually want to apply Eq. ( I36l) in order to get an expression of 
the shift of the front, we expand the latter formula for a, a <^ 1. The leading 
order reads 





(42) 




(44) 



^ [7r\2k-l)^ + 2\nq\ e 




(45) 



2 g(-lng)5/2 L"^ 



k=l 
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The shift of the front due to a fluctuation is obtained from ips with the help 
of Eq. ([36D: 



R{t,S) 



InU 



7o 



(7rx"(7o)t)^/2 
+00 

k=l 



{2k - 1) 



x"ilo)t 
L2 



2x"(70)« 



(46) 



q is the function of t given by Eq. fl39|) . This formula is extremely useful, since 
the series indexed by k converges fast. Even for moderately large values of t, 
a few terms accurately describe the whole function. This is actually the best 
formula for numerical evaluations of R. 

We shall now examine the limit of small t (y ^ 1). Then only the term k = 1 
has to be kept. The expression for R boils down to 



R{t, 5) = — In 1 + C2 
7o V 



(7rx"(7o))^/2 t5/2 



(47) 



As a final remark, let us note that the Poisson summation fl42p that we have 
used to rewrite the series of harmonics corresponds to a Jacobi identity for 
the d functions [24J. Equation f H6|) results from Eq. fHOj) with the replacement 

- = ^^pi^^ E (^'(2^ - 1)' + 2 Ing) e^^. (48) 



c?.^ Analytical expression for the correlations 



With the elements presented in the previous sections, we can write the ex- 
pression for cr^2 = ((-^1 — -^2)^)- It is enough to insert the expression for 



the probability of fluctuations (Eq. (fT6!) ) 
(Eq. (gUD) into Eq. (|2T|): 



and for the time-dependent shift 



dt 



2Ci 



I-C2 



e7o5 
2L3 



(49) 



where for 9gi?4(0|g) we use either one of the equivalent expressions ( HTl) . 
according to the limit that we want to investigate. We now have to fix the 
value of L. In Ref. [23], L was taken to be a constant. (The phenomenological 
model predicted L = Lq = lniV/70, but empirically, we saw that it was better 
to add a subdominant correction, namely L = Lq + + const.) In this 

case, a change of variable can be made in the integrand. All the parameters 
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may be factored out, leaving us with a simple numerical integral to perform: 







2C(2) = y. 



(50) 



Thus 



(51) 



dt 3jIL^ 

Replacing the product of the unknown constants by Eq. (|T8l) and q by Eq. ( I39l) 
and integrating over the time variable between and At = t — t/^h, we arrive 
at a parameter-free expression for function of At, namely 



2 
(^12 



271^ 



1 do 
2L^ q 



(52) 



We now investigate the two interesting limits, i.e. At ^ and At <^ L^. For 
large At, the integral is dominated by the region g — )■ 0, thus —dg'd^^Olq) may 
be replaced by its value at g = (— 9g't?4(0|0) = 2). Performing the remaining 
integration, we get 



^12 ~ — „ — At, (5-3) 

which is twice the second-order cumulant of the position of the FIP front, see 
Eq. ([111). For small At instead, say L ^ At ^ L^, we use the expansion of 
9gi?4(0|g) for g — )■ 1, i.e. the first term in Eq. fj48|) . which reads 

Wlg) = -V g(-lng)v' ^'^- ^54) 
Equation (l52l) boils down to the following expression: 



cr 
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27r3 / L2 \ 



At«L2 37o'\ X"(7o)At ^ V 2x"(7o)At^ 



So far, we have chosen the size of the front L constant, of the order of Lq. 
Another possible model for L would be to promote it to a function of S at the 
level of Eq. (149]) , namely 

L = Lq + 6 + const, (56) 

where the constant has to be determined empirically. This choice takes maybe 
into account more accurately the extension of the front by 6 that generates 
the fluctuations. The ^-integral cannot be performed analytically in Eq. (149!) 
except in some limits, so a priori, there is no simpler expression than Eq. (l49l) . 
Thus we need to know the values of Ci and C2 individually. We can consider 
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that Ci = 7o is the natural normahzation of the probabihty distribution p{S). 
Then, we must set C2 = 7r^x"(7o)/7o in order to satisfy Eq. (fT8l) . 

The above-mentioned two models, in which L is either constant or 5-dependent, 
differ by subleading terms in the large-L limit. Since the values of 6 which dom- 
inate the (5-integral in Eq. fH9|) are of order IuLq, like the first correction to 
Lq in the case of constant L, the models are not expected to differ significantly. 
We will check this statement numerically. 



3. 5 Scaling 



Looking back at Eq. ( 152|) . we see that has a nice scaling property. Indeed, 
we may rewrite the latter equation as 



2 

(^12 



D 



1o{vq - v) 



-Toc^o-")^* q 



(57) 



in terms of the properties of a single front (its velocity v and the diffusion 
constant D whose analytical expressions were given in Eq. ( !T2|) and f|T^ ). 
where vq can be read in Eq. (JTj). In particular, we have the following scaling: 



^12 



function[(t>o — f )At]. 



(5J 



From Eq. (!55|) . we see that the function in the right-hand side is exponentially 
damped when its argument is smaller than 1, i.e. parametrically for At -C L^. 

Once one knows the characteristics of the traveling waves in the FIP model 
(i.e. V and D), this scaling of the correlations is a pure prediction. Thus it will 
be interesting to check it in the numerical calculations. 



3. 6 Limits on the validity of the calculations 



Let us try and evaluate the limits on the validity of our calculations. The latter 
were essentially based on the assumption that the eigenvalue 7 = 70 of the 
kernel x dominates. While this statement is clearly true at large times, when 
the traveling- wave front is well formed, (see e.g. Ref. [21]), it must break down 
at early times right after a fluctuation has occurred: Indeed, a fluctuation has 
an initial shape that is far from the one of the asymptotic front, see Fig. [2J 

We wish to estimate the order of magnitude of the dispersion of the relevant 
eigenvalues about 70. To this aim, neglecting for the moment the boundary 
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conditions and the prefactors, we write the solution of Eq. (122|) as 

M(At,fc) ~ yd7e-^'=+^(^)^*. (59) 

The interesting values of k are the ones around the position of the wave front, 
therefore we write k = voAt + 6k, where 5k is oi the order of the size L of the 
front. Expanding xil) about 70, we write 

u{At, VoAt + 5k) ~ 6-^°^^^= J d{5-f) e-^7x5fc+ix"(7o)(57)^At+...^ ^gQ^ 

where ^7 = 7 — 70. It is clear from this equation that the relevant values of ^7 
are of the order of 5k/{x"{'lo)At). Since the order of magnitude of 5k is the 
size L of the front, we would a priori conclude that the dispersion of 7 around 
7o is small and hence that the calculation is valid as soon as At ^ L. 

However, we have also expanded xil) to second order. This means that for 
a generic kernel we have neglected terms of the form ^x^^K'lo)i5'y)^ At ~ 
L^/(At)^ (which would fit in the dots in Eq. (160|) ). The expansion is a good 
approximation if the latter term is small, i.e. if 

At > L^/^. (61) 



3. 7 Back to impact-parameter space 

So far, we have been working with the minimal model, consisting in two real- 
izations of the FIP model which evolve in the same way until their common 
tip reaches k^b, and which decorrelate for k > k/^h- The only relevant param- 
eter which determined the decorrelation of the positions of the fronts of the 
realizations was the time At = t — tAb after the tip had reached k^b- We now 
wish to discuss the transcription of the obtained results to impact-parameter 
space, which was our initial problem. 

To this aim, we will of course make use of Eq. ( !T5|) to express t^b with the 
help of the mean front velocity v. But we also need a length scale to which the 
distance in impact-parameter space Ab may be compared. The natural length 
is the dipole size at the position of the front, namely 

/,(t) = 2-^W=/,(tA6)2-^^*. (62) 

On the other hand, according to Eq. (fTSjl and disregarding the integer part 
operator, — logg Ab = k^b and the tip of the front k^b is ahead of the bulk 
^{tAb) by L: k^b = -^(^Afe) + L. Using the previous equation, we may now 
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express At as a function of Ab and of the length scale ls{t): 



At = - L + lOgr, 



Ab 



(63) 



V 



ls{t) 



The scaling fl58|) reads 




X function 




(64) 



This formula, together with the behavior of the scaling function (see Eq. flS^ ). 
shows that there is little 6-dependence until log2(A6/Zs(t)) ~ L^, that is to 
say, until Ab ~ lg[t)e^°'^^^^^^ . In other terms, the size Ab of the domain around 
impact parameter b in which the fluctuations in the position of the fronts are 
negligible is, in notations more familiar to QCD experts. 



where Qs{b) is the usual saturation momentum at impact parameter b. Note 
that since the fronts are statistically independent as soon as Ab x Qs[b) > 1, 
this result may seem a bit surprising: It says that the effective correlation 
length between different points in impact-parameter space is much larger than 
1 / Qs{b) in the parametrical limit of small a^. This is the main qualitative result 
of this paper. 



4 Numerical simulations 

In this section, we confront our analytical calculations to numerical simulations 
of the toy model. First, we consider the full model and test the validity of the 
assumption that the minimal model is a good approximation to the full model 
also for /3 ~ 1, i.e. when splittings to larger-size dipoles are authorized. Second, 
we compare the minimal model to the analytical results for the fluctuations 
between different positions in impact-parameter space (given essentially by 



4-1 Full model 

The model defined by Eqs. ([2]) and ([3]) is straightforward to implement nu- 
merically in the form of a Monte-Carlo event generator. The simplest is to 
store the number of dipoles in each bin in an array whose index i is related to 



e' 



,constxln^(l/a^) 



Ab 



Qsih) 



(65) 



Eqs. (in]),(E21)). 
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Fig. 3. One event of the full model with a = (3 = 1, N = 100 and t — t/\i, = 4. Only 
the bins A; > 5 are represented. (The bins for A; < 5 all contain N dipoles.) The 
number of dipoles in each bin is proportional to the blackness which is displayed. 
We see that in the transition region close to blackness, nearby bins are often of 
similar grey levels, which illustrates the statement that the density of gluons varies 
significantly only over scales which are larger than the relevant length scale ls{t) 
(see Eq. dM])). 

k and j through i = 2^~^ + j. The splitting dynamics relates bin i to 2i (down 
left), 2i + 1 (down right) and [i/2] (up; the square brackets stand once again 
for the integer part). 

We have to deal with an array whose size grows exponentially with time. It 
is thus very difficult to pick large values of t, and thus also large values of 
N. Indeed, the relevant time scale grows with N like In^ A^, and consequently 
the minimum number of entries in the array one wants to consider grows like 
gin N ^ practice, we limit ourselves to t < 4 and < 100. As for the time 
step dt, the most convenient is to take it small but finite. We set dt = 10~^. 

We start with one particle and evolve it for a few hundred units of time using 
the FIP version of the model. We obtain a traveling wave front, whose tip we 
eventually label k = 1. (The complete front sits in the bins A; < 1). From the 
initial condition built in this way, we evolve all bins for which A; < 1 using the 
FIP model, and all bins for A; > 1 using the full model. One event in shown in 
Fig. [3l Although A^ and t are small in this calculation, we see that the regions 
in impact-parameter space which have similar numbers of dipoles are larger 
than the local length scale ls{t) (see Eq. flB21) ). 

After the evolution times t = 3 and t = 4 respectively, we measure the position 
of the front at various impact parameters on a uniform tight grid ranging from 
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-| to +i. We use the following definition of the position of the front: 

fc=fco+l 

where ko is the largest /c for which ^^(A:,[Afex2'=-i]) = Note that in principle, 
we could have chosen X{Ab,t) = k^. In practice however, because of the 
discreteness of k in our model, this choice would introduce artefacts which we 
do not expect in real QCD. 

We compute the squared difference of the front positions between the impact 
parameters — | and — | + A6, and average over events. We plot the result as 
a function of t + log2 Ab/v, where v is the average front velocity measured at 
impact parameter — |. 

We compare the results to the correlations obtained in the minimal model, i.e. 
when we consider two independent realizations of an initial front. We do not 
attempt to compare to our analytical formulas since the values of N that we 
are able to reach are too small for the approximations that we had to assume 
to be relevant. 

The corresponding plot is displayed in Fig. H] for = 100, a = 1, /3 = 0, 
and in Fig. |5]with the same parameters except /3 = 2. First, we see that in 
the full model, the graph of cr^2 exhibits steps, i.e. (T^2 is constant by parts. 
This is related to the hierarchical structure of the model: The correlations 
between b = and any of the points at 6 > are identical; The same is 
true for —0.25 < b < 0, —0.375 < b < —0.25 etc... The logarithmic 6-scale on 
the t-axis makes the widths of the steps all equal. Next, we see that for small 
t — t/^b (i.e. impact parameters close to — |) there are very little fluctuations 
in the front positions. 

Finally, we see that for /3 = 0, as anticipated, the full model and the min- 
imal ones coincide almost perfectly (Fig. H]). For /3 = 2, i.e. when splittings 
towards larger dipole sizes are switched on and therefore new correlations ap- 
pear beyond the ones taken into account in the minimal model, there are some 
quantitative differences for large t (Fig. [5]). But we see that using the minimal 
model instead of the full model that keeps all impact parameters is a good 
approximation. This corroborates the conclusions of the work in Ref. [T8] . 

4-2 Minimal model 

We now set /3 = 0, in which discussed earlier and as checked numer- 

ically, the model exactly reduces to a collection of one- dimensional FKPP- 
like models. Hence, in order to compute two-point correlation functions, it is 



23 



TV = 100, (3 = 
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Fig. 4. £7^2 — ((^1 ~ ^2)^) as a function oi At = t — t^b in the full model with a = 1 
and /3 = (lines with steps; one corresponds to an evolution time t = 3, the other 
one to t = 4) and in the minimal model. In the full model, tAfe = [~log2^^]/^) 
where v is the measured velocity of the front at impact parameter — ^ . In the FIP 
model, t/\b is a fixed time, and corresponds to the time at which the tip of the front 
reaches k/\f,, the bin after which two uncorrelated evolutions take place. 
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Fig. 5. The same as in Fig. |3]but for /3 = 2. 
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N = IQi" 




Fie. 6. function of At = t — t/^t, in the minimal model with (3 = for 

A*" = 10^*^. We display the results obtained within the model in which the realizations 
decorrelate in the bins k > k^f, (labelled "Monte Carlo"), and within the model in 
which the decorrelation is complete after time tAfe (labelled "without correlations"). 
The theoretical curves use Eq. (I52p with the two possible choices for the front size 
L. Inset: The same, as a function of 1/At in order to highlight the small- At region 
where, as expected, important differences appear between the models. 

enough to evolve two realizations of the corresponding FIP model with the 
constraint that all bins with k < fc^fe be identical between the two realizations, 
and the bins k > k^f, be completely independent. Alternatively, we could also 
generate one single realization and evolve it for time steps, replicate it at 
time t^h, and then evolve the two replicas completely independently of each 
other. The difference between these two possible implementations of the min- 
imal model cannot be accounted for in our analytical calculations, thus the 
differences that we shall find numerically will give an indication of the model 
uncertainty. This time, our aim is essentially to check our analytical formulas, 
thus we will pick very large values of A^, even if they appear to be unphysical 
in the QCD context since they would correspond to exponentially small values 
of the strong coupling constant a^. 

The parameters of the model are obtained from Eq. ([8]) with a = 1, (3 = 
and dt = IQ-^: 

7o = 1.0136- ■■ , Wo = 2.6817--- , x"(7o) = 2.6098 ■■ ■ (67) 

These values are close to 1, e and e respectively, which would be the correct 
parameters if dt were infinitesimal, in which case ^(7) = (see Eq. ([H])). 
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Fig. 7. The same as in Fig. El for N = 10^^. All curves coincide almost perfectly. 

The numerical results are shown in Fig. El for = 10^*^ with the two versions 
of the model (we generated about 10^ realizations), and compared with the 
analytical predictions. We test the two possible choices for the size L of the 
front: Either L is a constant, which from our previous experience with FKPP 
traveling waves [23], we set to 

L = — lniV+ — InlnA^- — (68) 
7o 7o 7o 

or it is 5-dependent, namely 

L = —InN + 6 -lA. (69) 
7o 

The numerical constants, which are not determined in our theory, were chosen 
empirically so that they properly describe all numerical data for > 10^*^. 
In the first case, Eq. ( 152|) is used. In the second case, Eq. ( l49i) is integrated 
numerically over t and 6. We see that the agreement between the numerical 
calculation and the analytical predictions is good, except maybe for very small 
values of At where the calculations are not expected to be accurate. Indeed, for 
the same values of At, we also see in Fig. [6] a sizable discreapancy between the 
two versions of the minimal model. The calculations for = 10^° are shown 
in Fig. [71 The numerical results and the theoretical expectations (Eq. ( 1521) ) 
coincide almost perfectly. 

Finally, we check that the scaling in Eq. is well reproduced by the nu- 
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Fig. 8. Numerical check of the scahng ()58p . The curves for the different values of 
are very close together for > 10^'', but the scaling seems to break down for low 
values of N (see the curve for = 100), as expected. 

merical data. The Monte-Carlo simulations are shown in Fig. [HI plotted in the 
appropriate scaling variables. The diffusion constant of a single wave front D 
as well as the velocity v are measured from the same data. We see that all 
curves nicely superimpose for > 10^*^ (we show data for values of A^ as 
large as 10^°), while there are clear deviations for smaller A^ (see the curve for 
A^ = 100), as expected. 



5 Conclusion 



In this paper, we have built a model that possesses the main features of the 
QCD dipole model including the dynamics in impact-parameter space, and 
which is furthermore very easy to implement numerically. We have obtained 
analytical expressions for the fluctuations of the (logarithm of the) saturation 
scale from one position in impact-parameter space to another one nearby, 
which gives an indication on the homogeneity of the gluon number density. 

Since our analytical calculations are only based on some rather general prop- 
erties of the model, they should go over to full QCD after appropriate re- 
placement of the parameters, hopefully giving the correct small-a^ (large A^) 
asymptotics. 
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We have found that the saturation scale varies quite slowly. In the usual no- 
tations of QCD, if at position b in impact-parameter space the local satura- 
tion scale is Qs{b), then the saturation scale is uniform over a region of size 
gConstxin2(i/a2)/g^(^5^ arouud that position. 

Our calculations suffer the usual limitations in this kind of models: Analytical 
expressions are derived for exponentially large N, that is to say, exponentially 
small values of the strong coupling constant a.^- Although they often may be 
successfully extrapolated down to ~ 100 (a^ ~ 0.1), this is at the cost of 
tuning constants, and so far, we have not found a systematic procedure to 
compute finite- A^ corrections. 

The next step would probably be a numerical calculation in the case of full 
QCD, using, as proposed earlier [25], a combination of a Monte Carlo im- 
plementation of the color-dipole model in the low-density regime (such an 
implementation is already available; see Ref. [25|27|28] and Ref. f25] for more 
recent work) and of a numerical solution of the BK equation at the transition 
to saturation [17J. This looks very challenging, but maybe the nice smoothness 
of the saturation scale in impact-parameter space that we have found in this 
paper will help. 

Finally, it would be interesting to find an observable which would directly 
be sensitive to the gluon density at two points in impact-parameter space 
simultaneously. 
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